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Abstract. We describe a new algorithm for computing the ideal class group, 
the regulator and a system of fundamental units in number fields under the 
generalized Riemann hypothesis. We use sieving techniques adapted from the 
number field sieve algorithm to derive relations between elements of the ideal 
class group, and p-adic approximations to manage the loss of precision during 
the computation of units. This new algorithm is particularily efficient for 
number fields of small degree for which a speed-up of an order of magnitude 
is achieved with respect to the standard methods. 



1. Introduction 

Let K = Q{9) be a number field of degree n and discriminant A. In this paper, 
we present fast methods for computing the structure of the ideal class group of the 
maximal order Ok of K, along with the regulator and a system of fundamental 
units of Ok- 

Class group and unit group computation are two of the four principal tasks for 
computational algebraic number theory postulated by Zassenhaus (together with 
the computation of the ring of integers and the Galois group). In particular, they 
occur in the resolution of Diophantine equations. For example, the Pell equation 

T^-AU^ = 1, T,UeZ, 

boils down to finding the fundamental unit in a real quadratic number field of 
discriminant A (see [18 ). In addition, the Schaffer equation 

y2 = + 2*^ + . . . + (.T - 1)'^', fc > 2, 

can be solved using solutions to the Pell equation [T7]. Unit computations are key 
ingredients in solving almost all Diophantine equations, for example when solving 
Thue equations 5 . On the other hand, the computation of the ideal class group 
C\{Ok) of a number field K allows in particular to provide numerical evidence in 
favor of unproven conjectures such as the heuristics of Cohen and Lenstra [9] on the 
ideal class group of a quadratic number field, Littlewood's bounds [21] on L{l,x), 
or Bach's bound on the minimal bound B such that ideals of norm lower than B 
generate the ideal class group. The class group enters also into the computation of 
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the Mordell-Weil group of elliptic curves with the descent method, or the Brauer 
group computations for representation theory [TO] . 

In 1968, Shanks [HI [21] proposed an algorithm relying on the baby-step giant- 
step method to compute the structure of the class number and the regulator of a 
quadratic number field in time O (lAI^/"*"*"'), or O (|A|^/^+'^) under the extended 
Riemann hypothesis [20]. In 1985 Pohst and Zassenhaus [23^ published an algorithm 
that could determine the class group of arbitray number fields. Then, a subexpo- 
nential strategy for the computation of the group structure of the class group of an 
imaginary quadratic extension was described in 1989 by Hafner and McCurley [14j . 
The expected running time of this method is bounded by La %/2 + o(l)) where 

Li^{a,P) :=eW°g|Ar(iogiog|A|)i-"_ 

Buchmann generalized this result to the case of an arbitrary extension, the 
heuristic complexity being valid for fixed degree n and A tending to infinity. In a 
recent work |4], Biasse described an algorithm achieving the heuristic complexity 
certain classes of number fields where both the discriminant and 
the degree tend to inifinity. 

In parallel of theoretical improvements, considerable efforts have been invested 
to make the implementations of the subexponential methods efficient. In the qua- 
dratic case, Jacobson [16] described an algorithm based on the quadratic sieve for 
deriving relations between elements of G\{Ok)- He successfully used it for com- 
puting the class group and the fundamental unit of quadratic number fields. His 
implementation contained some of the practical improvements described in the con- 
text of factorization such as self initialization and the single large prime variant. 
This strategy was later improved by Biasse [3] who used a double large prime vari- 
ant and a dedicated Gaussian elimination technique. Attemps have been made 
to generalize sieving techniques to general number fields |22j . but the proposed 
algorithms remain impractical for the sizes of discriminant of interest. 

Our contribution. In this paper, we present an algorithm based on sieving tech- 
niques adapted from recent implementations of the number field sieve ^19j for com- 
puting C\{Ok) under the generalized Riemann hypothesis (GRH). We also describe 
a p-adic method for computing the regulator and a system of fundamental units. 
We show that these methods allow a significant improvement for number fields of 
low degree over the current state of the art. 

2. Generalities on number fields 

Let K he a number field of degree d. It has ri < d real embeddings {(Ji)i<r-i and 
2r2 complex embeddings {(Ji)ri<i<2r2 (coming as r2 pairs of conjugates). The field 
K is isomorphic to Ok ® Q where Ok denotes the ring of integers of K . We can 
embed X in /Cr if «) R ~ M''^ x C^, and extend the cr^'s to K^. Let T2 be the 
Hermitian form on i^R defined by 

T2{x,x') ai{x)a~i{x'), 



and let := ^JT2{x, x) be the corresponding L2-norm. Let {ai)i<d such that 
Ok = (Bi^Oii, then the discriminant of K is given by A = det^(T2(a,,aj)). The 
norm of an element x G A' is defined by N{x) = Hi ^^(s;)!. 
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To construct the ideal class group of Ok, we rely on a generalization of the 
notion of ideal, namely the fractional ideals of Ok ■ They can be defined as finitely 
generated Ojc-modules of K. When a fractional ideal is contained in Ok, we refer 
to it as an integral ideal, which is in fact an ideal of Ok- Otherwise, for every 
fractional ideal / of Ok, there exists r G Z>o such that rl is integral. Fractional 
ideals are assumed to be non-zero. The sum and product of two fractional ideals 
of Ok is given by 

IJ = {hji + + 1 1 e e l,ji,---ji e J} 

l + J = {i + j\tel,jeJ}. 

The fractional ideals of Ok are invertible, that is for every fractional ideal /, there 
exists := {x & K \ xl C Ok} such that 11^^ = Ok- The set of fractional 
ideals is equipped with a norm function defined by the index Af{I) '-= [Ok ■ I] 
for integral ideals, and the natural extension M{I / J) := J\f{I)/J\f{J). The norm 
of ideals is multiplicative and for principal ideals, it agrees with the norm of the 
generator N[xOk) = N{x). 

The ideal class group of Ok is defined by C\{Ok) '■— I-IV, where I denotes the 
group of fractional ideals of K and V Q I is the subgroup of principal fractional 
ideals. We denote by [a] the class of a fractional a in G\{Ok) and by h the cardinality 
of C\{Ok)- Elements of T admit a unique decomposition as a power product of 
prime ideals of Ok (with possibly negative exponents) . An element x g Ok is said 
to be a unit if {x)Ok = Ok, or equivalently if N{x) = 1. The units of Ok form a 
multiplicative group of the form 

U = nx (71) X ••• X (7^), 

where /x is the torsion subgroup of U, r := ri + r2 — I and the generators 7^ of 
the non-torsion part are called a system of fundamental units. The regulator is an 
invariant of K which allows us to certify the calculation of C\{Ok) and U. It is 
defined as i? = Vol(r) where F is the lattice generated by vectors of the form 

(ci log|7i|i, • • • ,Cr+ilog|7i|r+i), 

with \x\i := |cri(a:)| for i < r + I, ci = I for i < ri, Ci = 2 otherwise (for each 
complex embedding (Ji, if i < r + 1, then (fi = aj for some j > r -I- 1). 

3. The subexponential strategy 

The idea behind the algorithm of Buchmann is to find a set of ideals B = 
{pi, • • • ,Pn} whose classes generate C\{Ok), and then consider the surjective mor- 
phism 

— ^ / — ^ C\{Ok) 

(ei,...,ejv) > > njP.]^- 

From the fundamental theorem of algebra, the ideal class group satisfies C\{Ok) — 
Z''^/ker(7r o (p). Therefore, the knowledge of ker(7r o t^), which has the structure 
of a Z-lattice, enables us to derive C\{Ok)- In the meantime, elements of ker((/5) 
give us units as power-products of relations. From these units, we hope to derive 
a system of fundamental units of Ok ■ The subexponential strategy can be broken 
down to three essential tasks: collecting relations, calculating the class group and 
calculating the unit group. The subexponentiality is a consequence of a careful 
choice of B. 
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3.1. Relation collection. A preliminary step to tire relation collection is the 
choice of a generating set B = {pi, • • • ,Pn} of C\{Ok)- We choose the set of prime 
ideals of norm bounded by an integer B. The use of the Minkowski bound certifies 
the result unconditionally, but it causes the algorithm to take a time exponential 
in the size of A. To achieve subexponentiality, many authors chose the bound of 
Bach IJ, who proved that under GRH, C\{Ok) was generated by the classes of the 
prime ideals p satisfying Af{p) < 121og(|A|)^. Although asymptotically better, in 
practice this bound can be larger than the one described by Belabas et al. [5] who 
stated that under GRH, the class group was generated by the classes of the prime 
ideals of norm bounded by B satisfying 



^ logAA(p) 

(m,p):A/'(p'")<S ^ 



log7V(p'" 
log(S) 



> 2 log 1^1 



1.9n- 0.785ri 



2.468n+ 1.832ri 



log(B) 

In the rest of the paper, we assume that B is constructed with the bound of Be- 
labas et al. Indeed, Bach's bound enlarges the dimensions of the matrices that are 
processed during the computation of C\{Ok), thus inducing a slow-down that is 
not compensated by the fact that the relations are found more rapidly. 
During the relation collection phase, we collect relations of the form 

where (pi € K. We progressively build the matrix M :— (eij) G Z''^^ where k is 
the number of relations collected so far. Let A C ker(7ro(p) be the lattice generated 
by the rows of M. Operations on the rows of M allow us to retrieve a basis for 
A and its determinant. To determine if A has rank N, we perform operations 
modulo a random wordsize prime p. In particular, the LU decomposition of M 
modulo p allows us to identify the prime ideals that do not contribute to the rank 
of A. Additional relations involving these primes increase the rank of M, whose 
rows eventually generate a finite index sublattice of ker(7r). To find this index, we 
compute the Hermite normal form (HNF) of M, that is, we perform unimodular 
operations encoded by U G GLk (Z) such that 



/ 



UM 



hii 











V 



(0) 



with yj<i: < hij < hjj and Vj > i : hij — 0. Once the HNF of M is computed, 
adding new rows can be done very efficiently. In the meantime, the product J^,- hi^i 
gives us an indication on [A : ker(7r o (p)], as we see in § 13.31 



3.2. Class group computation. Given a matrix A G Z^^^ whose rows generate 
ker(7r o (^), unimodular transformations on both rows and columns of A yield the 
structure of CI{Ok)- More precisely. For every non-singular matrix A £ Z^^^, 
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there exist unimodular matrices U,V £ ll^^^ such that 

S := UAV = diag(di,...,dAr) 

where Vi such that 1 < i < : di+ijfii. The matrix S is cahed the Smith normal 
form (SNF) of A. 

Theorem 1. If the rows of A G 1/^^^ are a basis for keT{Tro(p)j and i/diag(di, . . . ,dpf) 
is the SNF of A, then 

C\{Ok) ^ Z/diZ X • • • X Z/dN- 

Once enough relations have been found, the rows of M generate ker(7r o ip)^ and 
the N non-zero rows of the HNF of M are a matrix A € Z^^^ whose rows are 
a basis for ker(7r o (p), and the SNF of A gives us C\{Ok)- However, finding the 
structure of C\{Ok) can also be done computing the SNF of a matrix which is in 
practice significantly smaller than A, namely the essential part of A. Indeed, for 
each matrix H in HNF, there exists an index I such that Vi > I, hi i = 1. The 
upper left I x I submatrix of H is called its essential part. As the classes of pi for 
i > I are generated by those of the pj, j < I, the SNF of the essential part of A 
suffices to recover C\{Ok)- 

3.3. Regulator and fundamental units computation. Computing the regula- 
tor and a system of fundamental units of K consists of finding kernel vectors of M. 
Indeed, if X = (xi, . . . , Xk) satisfies XM — 0, then we have 

In other words, 7 :~ Yii 4'V ^ unit. Every kernel vector X of M yields a unit, 
and we want to compute the group generated by all those elements as well as the 
regulator of this group, defined to be zero if the group is not of full rank. So far, 
finding of relations between units is mostly done using real linear algebra (LLL). 
The core problem here being the numerical instability of the matrices. This in 
itself is a consequence of the well-known fact that units are very large in general, 
writing the fundamental unit of a real quadratic fiels explicitly with the canonical 
basis needs exponentially many digits while it is always possible to find a product 
representation of size polynomial in log |A|. At the end of the procedure, we verify 
that the assumption we made on the completeness of the lattice of relations is true. 
To this end, we use an approximate of the Euler product 

where (k{s) — J2a a/ {a)" ^^'^ usual ^-function associated to K. Indeed, it allows 
us to derive a bound h* in polynomial time under ERH that satisfies h* < hR < 2h* . 
If the values det(r) and dct(A) do not satisfy this inequality, then we need to collect 
more relations. 

4. Sieving techniques 

In this section, we describe sieving techniques to derive relations in C\{Ok) 
for general number fields. This is a generalization of Jacobson's results [16] for 
quadratic number fields. We provide numerical data illustrating the considerable 
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impact of these techniques for class group and unit group computation in the case 
of low degree number fields. Given a generating set B — {pi, . . . Pn} for CI{Ok), 
the usual method for deriving relations consists of computing random exponents 
e := (d, . . . , Cat), a € Ok and a reduced ideal Ig such that 

pr---p%- = (a)/e. 

Then, every time is S-smooth, we obtain a relation. As the arithmetic of ideals is 
rather expensive when n > 2, the relation search in the computer algebra software 
Pari [28] and versions 2.x for a; < 18 of Magma [6] consists of enumerating short 
elements of Ig via the Fincke-Pohst method [12] . 

Our method consists of deriving relations from smooth values of polynomials, 
thus avoiding the cost of the ideal arithmetic and of the ideal reduction. Our method 
for finding smooth values is based on the recent development of the number field 
sieve algorithm [19]. The use of trivial methods such as trial division for finding 
smooth values of our polynomials would yield the same theoretical complexity, but 
would be impractical for large discriminants. The most efficient implementation of 
the enumeration-based strategy for finding relations is the one of Pari. Therefore, 
in the following, we assess the impact of our sieving method by comparing its 
performance with those of Pari. 

4.1. Polynomial selection. Let a be a ;B-smooth ideal of Ok- In this section, 
we show how to provide a polynomial Pa £ Z[X,Y] of degree n such that every 
{x,y) G 1? such that Pa{x,y) is i3-smooth yields a relation. Note that in theory, 
a can be any ideal, however, we obtained the best results by choosing a = Ok- 
Let a and /3 be two independent elements of o. Then, we create by interpolation a 
Pa e 1[X, Y] such that 

Va;,?/eZ2, Pa{x,y) =Af{xa + yl3). 

Every time (t)x,y ■= xa + yf3 has a smooth norm, we add the relation corresponding 
to the principal ideal {4>x.y) to the relation matrix. Before applying sieving algo- 
rithms to Pa to derive relations, we need to ensure that it is likely to yield enough 
smooth values. Polynomial selection is an important part of the number field sieve 
algorithm, and so it is in our algorithm. However, the specificities of our context 
prevent us from directly adapting the methods of NFS for selecting the sieving poly- 
nomial. First of all, we can afford to find relations with many different choices of 
a and /?, whereas the choice of a sieving polynomial in the NFS algorithm is fixed. 
We require that our choices of a and (3 yield polynomials with small coefhcients, 
and that we have a sufficient randomization at the infinite places to avoid drawing 
<f>x,y spanning the same subgroup of the unit group of Ok- 

To randomize the choice a, /3, we consider random coefRcients ai, . . . , a„ € M" 
such that X]i<n '^i ~ ^- -'^'-'^ every such n-tuple a, we define the embedding 

a — > M" 
il^S- a I — > (ailog|Q;|i,...,a„log|Q;|„). 

For every choice of a, the set of elements of the form ■(/'a(a) for a G a is a lattice 
As of M" for which we can find an LLL reduced basis for the norm 

2 ■■ [Xi,. . . ,Xn) ^ e ^x^ + ... + e "Xn. 

For every choice of a, the first two vectors a, (3 of an LLL reduced basis of As 
are potential candidates for the creation of a polynomial yielding smooth values. 
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Every time we draw such a couple of elements of a, we need to make sure that they 
do not generate the same Z-module as another pair previously used. To prevent 
this from happening, every time we draw a couple a, /3 by the previous method, 
we express them in terms of the canonical Z-basis of Ok- Thus, to every couple 
a, P corresponds the matrix Ma^p G 1?^" of their coordinates. The HNF of Mq,^^ 
uniquely represents the Z-module spanned by {a, 13). Thus, to avoid duplicates, 
we store a hash of the HNF of Ma^p in a hash table every time we use a couple 
(a, /3) to draw relations. We summarize the procedure of the selection of a sieving 
polynomial in Algorithm [1] 



Algorithm 1 Polynomial selection 
Input: a, {Ai, . . . , An) , HashTable 

Output: Sieving polynomial Pa,i3 corresponding to e o 
1: while a new a, f3 has not been found do 

2: Draw ai < Ai, . . . , a„ < A„ at random such that ai + . . . + a„ = 
3: Let a, (3 be the first two elements of a LLL-reduced basis of A^ for a — 
(ai, . . . ,an) 

4: Compute the hash ha^p of the HNF of Ma,p 
5: if ha^p ^ HashTable then 

6: Compute by interpolation Pa p ^ "^[^^ Y\ such that Pa p{x, y) — N{xa + 

vP) 

7: end if 

8: end while 

9: return a,l3,Pa.p 



4.2. Line sieving. The quadratic sieve algorithm [25] used to derive smooth values 
of a binary quadratic form generalizes to the case of polynomials of arbitrary degree. 
Its design follows from the observation that if P € Z[A, Y] is a polynomial of degree 
n, then 

(1) Vj/o e Z, p I P{rp, yo) ^ Vi G Z, p\ P{rp + ip, yo). 

Given yo G Z, we wish to find the x G [— //2,//2] such that P{x,yo) is i?-smooth, 
where B is the bound on the norm of the prime ideals in the factor base. Instead of 
trying them all, we prefer to isolate a short list of good candidates that we test by 
trial division, lip \ P{x, yo) for many p < B, then P{x, yo) is likely to be i?-smooth. 
From ([T|), we know that once we have one root Vp of P{X,yo) modp, then we can 
derive all the others by translation by {p,0). Line sieving consists of initiating to 
zero an array S of length / whose cells represent the x G [— //2,//2]. Then, for 
each p < B,we compute the smallest roots Xp G [—1/2, 1/2] of P{X, yo) mod p and 
repeat 

S[xp] ^ S[xp] + log(p), Xp ^ Xp +p. 

Then, whenever S[x] w log(P(x, yo)) for x G [—1/2, 1/2], the value P{x, yo) is likely 
to be P-smooth. We summarize this procedure in Algorithmic] 

4.3. Lattice sieving. Let Pa,i3{X,Y) G Z[X,Y] be the sieving polynomial de- 
scribed in § 14.11 B the bound on the norm of the ideals in the factor base, and 
I,J € Z>o. Every couple {x,y) G [— 7/2, //2[x [1, J] such that Pa,p{x,y) is B- 
smooth yields a relations. Therefore, one can repeat the line sieving operation on 
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Algorithm 2 Line sieving 

Input: P G Z[X, y], /, B,yoeZ 

Output: Smooth values of P(X,?/o) in [-1/2,1/2] 

1: L ^ 0, Vx e [-//2, 1/2] , S'ia;] ^ 0. 

2: for p < B do 

3: Let Xp be the smallest root of P{X, j/o)niod in [—1/2, 1/2]. 
4: while Tp < //2 do 

5: S[Xp] S'[Xp] + log(p), Xp ^ Xp+p. 

6: end while 
7: end for 

8: for X e [-//2,//2] do 

9: if S[x] w P{x,yo) then 
10: If P(a:;, yo) is B-smooth, L <— L U {a;}. 
11: end if 
12: end for 
13: return L 



Pa.piX, yo) for every yo e [!,</]• This method is efficient when sieving with primes 
p < I. but when the primes are significantly larger than /, the root computation at 
Step 3 of Algorithm [5] is often performed for nothing since there is a good chance 
that none of the x G [—7/2, //2[ will be a root of Pa,i3{X,yo) mod p. A way around 
that is to have an array S of length IJ representing [—1/2, 1/2[ and to fill it by line 
sieving methods for the primes p < I and by lattice sieving for the other primes. 

The lattice sieve was first described by Pollard [24]. Since then, it has been 
extensively studied and improved in the past 15 years, and the most recent devel- 
opments of this methods yeld the factorization of RSA768 (see [19]). This strategy 
relies on a one-time enumeration of roots of Pa.p{X, Y) mod p in [—1/2, 1/2[x [1, J]. 
The entry x < IJ of the array S that we use to store the logarithmic contributions 
corresponds to the couple £ [— 7/2, //2[x [1, J] where 

i = {x- 1/2) mod / 
]^{x-i- 1/2)/ 1. 

As in the line sieving case, every entry of S is initialized to zero, and for every p < B 
and every G [— //2, //2[x [1, J] such that p | Pa^i3{i,j), we want to perform 

the operation S[x] ^ S[x] + log p. Line sieving repeated on every line j < J allows 
us to efficiently do this for p < I. For the others, we followed the approach of |13j. 
as it is done in [19! for the factorization of RSA768. By [13j Prop. 1], we know 
that for every p such that we have a root of Pa^p{X, 1) modulo p, there exists a 
basis {(a, &), (c, d)} of the couples {i,j) such that p | Pa,p{i,j) that satisfies 

• & > and d > 

• -I <a<Q<c< I 

• c — a > I. 

This basis is computed via an algorithm described in |13^ that relies on the continued 
fraction expansion of Vp. To fill the array S, we start from {i,j) — (0,0) which is 
a common root modulo all primes. Then, by induction, we construct the next pair 
{i',j') from {i,j) by choosing 

• (*! j) + (flj ^) if « > —a 
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• {i, j) + {c, d) ii i < I — c 

• ih j) + ia^b) + {c, d) if I — c < i < —a. 

4.4. Special-^/. The sieving space //2, //2[x [1, J] only contains a limited num- 
ber of couples yielding a smooth value. Enlarging / and J might cause its size 
to rapidly exceed the single precision. Let g be a prime, the special-g strategy con- 
sists of sieving with a polynomial Pq derived from the original sieving polynomial 
P such that 

V(*,j) e [-//2,//2[x[l, J], 3{x,y) e I^P^.j) = Pi^.v) 
V(i,j)e [-//2,//2[x[l,J], q\Pq{i,j). 

This strategy was used by Pollard in his original paper [24] to sieve on the rational 
side, but most current implementations use it in the rational side as well [19]. 
To create Pq for a given q, we need a root of P modulo q. Then, we find a 
reduced basis (oq, 6o)) (oi; ^i) of the lattice spanned by the vectors (r^, 0), (r^, 1). 
The polynomial Pq is then simply given by 

Pqihj) = P{iaQ+jai,ibQ+jbi). 

The reduced basis is given by successive Gaussian reductions, as explained in [13] . 
Then, to sieve with a given polynomial P, we repeat the procedure described in ^4.31 
for many different polynomials of the form Pq. Fortunately, once the roots of P 
mod p for all p < B have been computed, it is possible to use these values to 
compute the roots of Pq mod p for p < B. Indeed, 

P{iao + jaijibo + jbi) = mod p 

means that there is some root Vp of P{X, 1) mod p such that Vp = mod p. 

This implies that we have Pq(r'^, 1) = mod p for 

i ai- rpbi 

= - = — mod p, 

J ao - rpbo 

which gives us a root of Pq{X, 1) mod p from (ao, bo), (ai, &i) and a root of P{X, 1) 
mod p. We summarize our procedure to derive relations from an ideal a C Ok in 
Algorithm [3l 

4.5. Overall relation collection phase. A necessary condition to compute the 
class group and the unit group is to produce a full-rank relation matrix M . Our 
sieving methods allow us to derive relations in CI{Ok) very rapidly, but it is hard 
to force a given prime to occur in a relation. The best performance is obtained 
by sieving with the trivial ideal Ok ■ If we want to see a given prime ideal p | (p) 
occur in a relation, one can use the special-g with q = p, or sieve with the ideal p. 
However, even after using those methods, some prime ideals still do not contribute 
to the rank of M. Rather than sieving in random power-products involving missing 
primes, one might prefer to switch to enumeration-based methods to complete the 
relation search. To identify the primes that need to appear in a relation, we perform 
an LU decomposition of the relation matrix modulo a random wordsize prime. We 
try to produce enough relations with sieving so that the rank of M is 97% of #K. 
Then we find additional relations with enumeration. 

To assess the advantage of sieving over enumeration techniques, we need to isolate 
its contribution to the performances of the class group and unit group computation. 
To do this, we used a modified version of the function bnf init of the computer 
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Algorithm 3 Sieving procedure 

Input: a C Ok, B = {p \ M{p) <B} ,I,J e Z>o. 
1: Select a, (3 G Ok and a sieving polynomial Pa.p with Algorithm [TJ 
2: Vp < B, compute the roots of Pa^piX, 1) mod p. 
3: for q < B do 

4: Compute Pq and its roots modulo the p < _B as in MAi 
5: Let S be an array of size IJ initialized to 0. 
6: for p < / do 

7: Do S[x] ^ S[x] + log(p) for each x representing G [— //2, //2[x [1, J] 

such that p I Pq{i,j) by repeating Algorithm [2] for each line j < J. 
8: end for 
9: for p > I do 

10: Calculate a basis {(a, b), (c, d)} of the lattice of points in [—1/2, 1/2[x [1, J] 
that are roots of Pq{X, Y) mod p with the method of H4.'S[ 

11: Do S[x] ^ S[x] + log(p) for each x representing (i, j) e [-1/2, 1/2[x [1, J] 
such that p I Pq{i,j) by using the method of M.3\ 

12: end for 

13: end for 

14: for a; < / J do 

15: if S[x] sa \og{Pq{i,j)) where x represent G [-//2, //2[x [1, J] then 
16: If log{Pq{i, j)) is _B-smooth, store the corresponding relation. 
17: end if 
18: end for 



Algorithm 4 Full rank relation matrix computation 
Input: K, B 

Output: A full-rank relation matrix for the primes of norm bounded by B 
1: B^{p\M{p)<B}^{pi,--- ,Pn}. 

2: Derive N relations by repeating Algorithm |3] with a = (1). Let M be the 
relation matrix. 

3: Perform an LU decomposition of M and let EmptyList be the list of zero 

columns. 
4: for p e EmptyList do 
5: Sieve with p, update M. 
6: end for 

7: Update EmptyList by updating the LU decomposition of M. 
8: for p G EmptyList do 

9: Find a relation involving p by enumerating short elements in random power- 
products. 
10: end for 
11: return M 



algebra software Pari that accepts in input a list of precomputed relations. We 
interfaced via SAGE this version of Pari with a developping version of Magma 
containing a function creating relations with the sieving algorithm. The Magma 
function tries to create enough relations so that the rank of M be 97% of 
and passes it to Pari which adds new relations with enumeration methods and 
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Table 1. Impact of sieving on class group and unit group com- 
putation of small degree number fields 



n 


loga |A| 


Pari 


Pari+Sieving 


3 


120 


76 


11 


3 


140 


694 


66 


3 


160 


6828 


333 


3 


180 


29807 


2453 


4 


120 


38 


7 


4 


140 


366 


24 


4 


160 


4266 


175 


4 


180 


31661 


1201 


5 


120 


33 


18 


5 


140 


295 


64 


5 


160 


3402 


378 


5 


180 


16048 


2342 


6 


120 


40 


111 


6 


140 


294 


161 


6 


160 


1709 


1012 


6 


180 


14549 


8413 



calculates the class group and the unit group. We compared the performance of 
this approach to the traditional bnf init function of Pari. There are two main 
reasons for using a hybrid version. The first one is that Pari's implementation of 
enumeration techniques is the most efficient. As these are necessary to finish the 
creation of the relation matrix after calling the sieving algorithm, it is interesting 
to see how the two perform together. Another reason for this choice is the fact that 
many different algorithms contribute to the computation of the class group and the 
unit group. In particular, we use time-consuming linear algebra methods such as 
the HNF computation. Our methodology avoids the risk of seeing the influence of 
the quality of the implementation of other algorithms occuring in the class group 
and unit group computation. 

We performed our computations on a 2.6 GHz Opteron with 4GB of memory. 
We used a branch of the developping version 2.6.0 of Pari provided by Loi'c Grenier 
and the developping version of Magma, interfaced via Sage 4.7.2. We allocated 3GB 
of memory to the computation made with Pari. For each size d, we drew at random 
10 number fields with discriminant satisfying log2 |A| = d. For each discriminant, 
we computed the class group and the unit group with bnf init, which we refer 
to as the "Pari" method, and with the hybrid version which we refer to as the 
"sieving-|-Pari" method. The average timings, in CPU sec (rounded to the nearest 
integer), are presented in Table [T] They illustrate the impact of sieving methods 
for small degree number fields. It is very strong for degree 3,4 and 5 number field, 
for which we often witness a speed-up of a factor at least 10 while it is rather 
moderate for degree 6 number fields, and negligible for number fields of degree 
7,8. Finding smooth values of a polynomial gets more difficult when we increase 
its degree, but it is not the only reason why the impact of sieving decreases with 
the degree. Indeed, for degree 6 number field, our sieving algorithm still derives 
relations at a competitive pace, but there are many linear dependencies whereas 
enumeration allows a more targeted search, thus avoiding linear dependecies. To put 
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Table 2. Impact of the quadratic sieve on fields generated by 
X2 + 4(10" + 1) 



n 


20 


30 


40 


45 


50 


54 


Magma 2.18 


0.7 


6 


22 


128 


170 


1453 


Pari + Sieving 


0.5 


5 


44 


271 


593 


1085 


Pari 


0.2 


3 


66 


556 


1562 


9251 



these improvements into perspective, we show in Table [2] the impact of Jacobson's 
self initializing quadratic sieve |16) which is implemented in Magma 2.18. The 
timings for "Pari" and "Pari + Sieving" are derived under the same setting as for 
Table [H In addition, we added the performances of of Magma 2.18 which uses 
different methods for linear algebra. Timings for the same serie of number fields 
were reported by Jacobson in [ini Table A. 3] on a 296Mhz SUN processor (for a 
fair comparison one has to take into account the verification time since the timings 
of Table [1] and Table [2] correspond to a certification under GRH) . 

5. Computing the unit group 

Assume that we have created a relation matrix (e^j ) corresponding to the rela- 
tions 

(0,)=pr'^---pr- 

Every kernel vector allows us to derive a unit of Ok- Let Pi^ - ■ ■ , /3fe be a generating 
system of the unit created so far. We compute a new unit /3', and we wish to find a 
new minimal generating set for (/3i, • • • ,/3fc,/3). Usually this is done by computing 
(real) logarithms of the units followed by some approximate linear algebra to find 
a (tentative) relation as well as the (tentative) new basis. This then is followed by 
an exact verification of the relation to guarantee correctness. The difficulty comes 
from the fact that the entries in the real matrix differ vastly in size - by several 
orders of magnitude - thus making it neccessary to work with a huge precision, in 
fact the precision is also subexponential in the discriminant for guaranteed results. 

Here, we propose to use p-adic logarithms instead. The key advantage comes 
from the much better control of error propagation in the linear algebra: unless 
division by non-units happens, linear algebra does not increase errors. However, 
while the correctness is based on the unproven Leopold conjecture about the non- 
vanishing of the p-adic regulator, this is not a problem in practice: any relation 
found by the p-adic method can easily be verified unconditionally, thus a failure of 
the algorithm would provide a counter example to Leopold's conjecture. 

We start by choosing a prime p such that the p-adic splitting field Kp has mod- 
erate degree, here we allow at most degree 2. Then we have n embeddings (pi 
oi K ^ Kp, and we define a map Lp : K* — > Kp : x H> {log (j)i{x))i where 

is the usual p-adic logarithm extended to Kp. In order to estimate the nec- 
essary p-adic precision, we also need the usual real logarithmic embedding, de- 
noted by L : K* — R''"'"^. We are looking for a (rational) solution (xi) to 
J2xiLp{/3i) = Lp{/3'). Using p-adic linear algebra we will instead get a p-adic so- 
lution (or a proof that /3' is independent). Using standard rational reconstruction 
techniques, we derive the rational solution from the p-adic one and then the inte- 
gral relation between the units. In order to estimate the p-adic precision, we bound 
numerator and denominator using Cramers rule and universal lower bounds on the 
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logarithms of units. The rational solution then also satisfies J^^i^iPi) — 
Let {ai)i be a basis for . . . , (3s, (3'). By Cramer's rule, we write 

X, = det(Lp(/3i), . . . , . . . , det(Lp(/3i), . . . , 

Since the (unknown) (a^) form a basis, we see that 

det(ip(/3i), . . . , Lp(/3'), . . . , det(ip(ai), . . . , 

is an integer and the same is true for L instead of Lp, thus we can write Xi as a 
quotient of integers. In either case, to make sense of the determinants, we will have 
to select an apropriate number of rows to make the matrices square. To bound the 
integers, we make use of the Hadraniat bound for det(L(/3i), . . . , £(/?'), . . . , i(/3s)) 
and some universal lower bound for det(_L(ai))i. For the lower bound we use lower 
bounds of logarithms of non-torsion units (jTTJ 3.5]): ||L(ai)||2 > or, if the 

unit group has full rank, s = r = 7'i + 7'2 — 1, we use lower regulator bounds, possibly 
coming from the Euler product. Having obtained bounds from the real logarithm 
(L) with low precision, we calculate the p-adic precision required to find Xi using 
p-adic linear algebra and rational reconstruction. In the course of the computation 
it can happen that the p-adic determinants (p-adic regulators) have non-trivial 
valuation. In this case we have to restart the computation with a correspondingly 
higher precision to account for the loss. Since the Leopold-conjecture is non-proven 
as of now, we also need to verify the solution by computing a low-precision estimate 
for II '^XiL{f5i) — i(/3')ll to compare it to the lower bound used above. 

From the relation Xi we can easily obtain a presentation of the new basis in 
terms of the /3i, /?'. For optimization, we then proceed to compute a new basis 
such that the real logarithms are (roughly) LLi-reduced. We note that we do not 
rely on any LLL estimates here, so any heuristic algorithm that aims at reducing 
the aparent size will do. Since we do not have any LLL algorithm that will accept 
real input (as opposed to rational), it is important that this does not influence the 
correctness. 

5.1. Advantages of the p-adic method. There are two core advantages of the p- 
adic logarithms over the ordinary, complex, ones: first, the linear algebra problems 
we need to solve in order to find dependencies or relations between units have a much 
simpler error analysis. In fact, contrary to the complex case, it is possible through 
the use of ring based operations to solve linear equations without any additional 
loss of precision. This is very important in the context of unit computation since 
the matrices representing the image of L{a) are very badly conditioned for classical 
numerical methods. The other advantage of the p-adic logarithms is more subtle: 
if we assume Leopold's conjecture to hold for the field(s) we are interested in, then 
instead of doing linear algebra over R with a precision of say q to find dependencies, 
it is sufficient to work with a real precision of q/2 and a p-adic precision of q/2 as 
well. Thus, assuming classical multiplication, we gain a factor of about 4 through 
the use of lower precision. Using fast multiplication (in high precision), the gain 
is smaller but still noticeable. But the most important advantage is the much 
easier precision control: instead of complicated and very delicate estimates for 
linear algebra problems, all we need are upper bounds on linear combinations with 
integral coefficients - which are trivial to obtain. 

We should also mention that one disadvantage of the p-adic method lies in the 
total lack of control over the real size of the units, thus it needs to be paired 
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with a crude (and uncritical for correctness) size reduction algorithm. Also, it is 
(currently) not possible to forego completely the use of complex (or real) logarithms, 
as the p-adic method is not capable to proving a unit to be torsion - without the 
knowledge of bounds on the real size. 

5.2. Lower bound from Euler Product. Suppose that, as in the class group 
algorithm, we are given an approximation of the Euler product, ie. we have a real 
number E such that l/\/2 < hR/E < \f2. After the relation matrix has full rank, 
and assuming the factorbasis is large enough for correctness, we have an upper 
bound for the class number, thus a lower bound for R. This lower bound will be 
several orders of magnitude larger than the universal bounds available otherwise. 

5.3. Saturation. After the inital steps of the algorithm, when the relation matrix 
has full rank, we have a tentative class number h and a tentative regulator R. 
Experimentally, at this point, liR does not approximate the Euler product very 
well - the product will be off by several orders of magnitude. However, after we 
found one or two more relations, the product has the same size than the Euler 
product, it frequently even looks like only a factor of 2 is missing in either h or R. 
To find the last missing relation can easily take more time than the entire previous 
run therefore we suggest using saturation methods instead. At this point in the 
algorithm the relations define a subgroup V of the S'-unit group Us where S is the 
factor basis. From the Euler product we know that the index [Ug : U) b is small, 
lets say b < B. For any prime p\b there is some u £ Us\U such that G U. Let 
us fix the prime p. For any prime ideal Q ^ S such that p\Af{Q) — 1 we can define 
the map (f)Q : U ^ ^*q/(^*qY mapping 5- units into the multiplicative group of the 
residue class field modulo p-th powers. The Cebotaref theorem |29] guarantees that 
if M € [/ is not a p-th power, there will be some Q such that 4>q{u) is non-trivial, 
ie. u is not a p-th power modulo Q. We now simply intersect ker(/)Q for several 
Q until either the intersection is or it is does not change for five consecutive 
Q. We expect that any u G U/ Ci kerc^g will have a p-th root in Us but not in S. 
Therefore = u is a new relation that will change hR by p. Repeating this for all 
p < B until we cannot enlarge U any more we find the missing relations. 

5.4. Representation. During the execution of the algorithm, all (S'-)units are 
naturally represented as power products of the relations coming from the sieving 
(or the saturation) . It is well known that the explicit representation of the units with 
respect to a fixed basis for the field can require exponentially large coefficients, so it 
is important to operate on the power products as much as possible. However, even 
the exponent vectors constructed for the basis of the unit group, or the saturation, 
will become huge, so we need to "size reduce" the power products. In particular, 
this happens even if the resulting element is not too large. Using ideas of ^ for 
compact representations and [15 for reduced divisors in function fields, we can find 
a representation for those elements that depends only on the logarithmic size (and 
the number field) rather than the execution path. For any prime p we can write 
any unit 

with elements such that the size of ai depends on the discriminant and p only. The 
length of the product comes from L{u). Furthermore, in this presentation it is easy 
to test for p-th powers as only ao needs to be tested and this is a small element. 



NEW TECHNIQUES FOR NUMBER FIELD COMPUTATIONS 



15 



5.5. Example. To illustrate the power of the p-adic method, we look at a totally 
real quartic field generated by a root of 

+ 17211a;^ + 5213a;^ - 176910463a; - 4958. 

The discriminant A of the maximal order has 38 digits. In the course of the 
computation, we found 534 relations involving prime ideals of norm up to 3000 = 
0.4 log'^ I A| describing a trivial class group. We then searched for 5 further relations 
to obtain units Ui {1 < i < 5). As power products of the relations, the units are 
given via exponent vectors with ||ei||oo ranging between 10®° and 10^^° and 
20 < ||ei||i/||ei||oo < 92. So, while not uniformely large, the exponents are non- 
sparse, involving huge integers. Using a decimal precision of 170 digits, wc establish 
that the logarithms of the units are roughly ||L('Ui)||oo « 10^^°. The first three 
units are indeed independent, giving a basis for a subgroup of full rank, the fourth 
is then dependent. Choosing the prime p = 10337 we get Qp as a splitting field. 
Using a p-adic precision of 245 digits (ie. working in Zp mod p^^^), we compute the 
dependency for the fourth unit, involving exponents of around 10^^°. The new unit 
group is then tentatively LLL reduced, producing a new basis where the 11^(^^)1100 
are bounded by 10^ only. The last unit then involves a much smaller dependency, 
here the exponents are only around 10^°. 

Unfortunately, looking at the Eulcr product, the unit group is not complete. 
However, the saturation technique outlined above takes Isec to determine that 
the product of the three basis elements is (proabably) a square. Finding a better 
representation where the exponents are all powers of 2 takes less than 1 sec and 
then we can enlarge the unit group easily. 

Due to the implementation, the p-adic precision used was actually higher: chang- 
ing (increasing) precision is very compuationally expensive, so we try to avoid this 
and simply double the precision. We used a precision of 320 for the p-adics and a 
maximal precision of 1000 for the real precision. The computation of the log is the 
dominating part, we spent 50sec or 90% of the total processing time here. 

6. Conclusion 

We introduced new techniques to enhance the performances of the subexponen- 
tial methods for computing the class group and the unit group of a number field. 
In particular, sieving allows a speed-up of an order of magnitude for number fields 
of small degree. These techniques could be developped even further. Indeed, we 
have not taken into account all the improvements to sieving techniques describes in 
the context of the number fields sieve algorithm such as large prime variations or 
cache- friendly methods. It is also notable that fast techniques for deriving relations 
in the class group of a small degree number field have applications in evaluating 
isogenics between small genus curves via complex multiplication methods. Indeed, 
in that case, evaluating isogenics between genus g curves involves relations in the 
class group of a degree 2g number field. 
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